Local Approximations to the Gravitational Collapse of Cold 

Matter 



Lam Hui and Edmund Bertschinger 
Department of Physics, Massachusetts Institute of Technology 
Cambridge, MA 02139 

ABSTRACT 

We investigate three different local approximations for nonlinear gravitational 
instability in the framework of cosmological Lagrangian fluid dynamics of cold 
dust. By local we mean that the evolution is described by a set of ordinary 
differential equations in time for each mass element with no coupling to other 
mass elements aside from those implied by the initial conditions. We first show 
that the Zel'dovich approximation (ZA) can be cast in this form. Next we 
consider extensions involving the evolution of the Newtonian tidal tensor. We 
show that two approximations can be found that are exact for plane-parallel 
and spherical perturbations. The first one ("non-magnetic" approximation, or 
NMA) neglects the Newtonian counterpart of the magnetic part of the Weyl 
tensor in the fluid frame and was investigated previously by Bertschinger & 
Jain (1994). A new approximation ("local tidal", or LTA) involves neglecting 
still more terms in the tidal evolution equation. It is motivated by the analytic 
demonstration that it is exact for any perturbations whose gravitational and 
velocity equipotentials have the same constant shape with time. Thus, the LTA 
is exact for spherical, cylindrical, and plane-parallel perturbations. We tested all 
three local approximations in the case of the collapse of a homogeneous triaxial 
ellipsoid, for which an exact solution exists for an ellipsoid embedded in empty 
space and an excellent approximation is known in the cosmological context. 
We find that the LTA is significantly more accurate in general than the ZA 
and the NMA. Like the ZA, but unlike the NMA, the LTA generically lead to 
pancake collapse. For a randomly chosen mass element in an Einstein-de Sitter 
universe, assuming a Gaussian random field of initial density fluctuations, the 
LTA predicts that at least 78% of initially underdense regions collapse owing to 
nonlinear effects of shear and tides. 

Subject headings: cosmology: theory — dark matter — gravitation — large-scale 
structure of universe 
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1. Introduction 

The complexity of nonlinear gravitational instability challenges our understanding of 
the universe. Even though the law of gravity between two bodies is very simple in the 
non-relativistic limit, the long-range interactions among exceedingly many bodies leads to 
behavior that defies simple analysis beyond the linear regime. Computer simulation with 
N-body methods provides a comprehensive approach to this problem, but it suffers from 
finite dynamic range and computational expense. Even more importantly, simulations do 
not increase our understanding of dynamics without guidance from analytical approaches. 

In this paper, we explore a class of what we call local approximations for the nonlinear 
dynamics of self-gravitating cold matter. By local we mean that each mass element behaves 
as if it evolves independently of all the others once the initial conditions are specified. This 
might sound quite implausible. After all, mass elements do influence each other through 
gravity. However, as we will demonstrate, the celebrated Zel'dovich (1970) approximation 
(henceforth ZA) can be viewed as exactly an approximation of this sort. 

In the past several years, there have been various attempts to improve upon the 
ZA, including the adhesion approximation ( [Kofman, Pogosyan, fc Shandarin 1990[ ), the 



frozen flow approximation flMatarrese et al. 1992j ), the frozen potential approximation 
( Brainerd, Scherrer, fc Villumsen 1993| ; |Bagla fc Padmanabhan 1994|) , the truncated 



Zel'dovich approximation QColes, Melott, fc Shandarin 19931) , an d higher-order Lagrangian 



perturbation theory ( [Melott, Buchcrt, fc Weiss 1995|) (note that the ZA can also be 
regarded as the first-order solution in Lagrangian perturbation theory). Most of them are 
attempts to deal with the evolution of high density regions after trajectories cross, when 
the ZA ceases to be adequate. However, this is a difficult problem. Aside from the spherical 
model ( [Peebles 1980| ) and its cousins, there still exists little in the way of approximations 
methods for post-collapse evolution. 

In this paper, we will not try to tackle the problem of trajectory crossing or the 
subsequent nonlinear evolution. Instead we ask whether one can improve upon the ZA 
even before orbits cross by seeking generalizations of the ZA within the framework of local 
approximations. In simple terms, a local approximation is one in which the evolution 
of each mass element is described by a set of ordinary differential equations in time in 
which there is no coupling to other mass elements, aside from those implied by the initial 
conditions. For instance, as we will explain more fully later, the evolution of a given mass 
element under the ZA is completely determined once the initial expansion, vorticity, shear 
and density at this mass element are specified. (The first three quantities correspond to 
the trace, antisymmetric part and traceless symmetric part of the velocity gradient tensor.) 
The evolution of other mass elements have no effect on the evolution of these quantities 
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at this mass element. In other words, under the ZA, all the information about other mass 
elements is encoded in the initial conditions. Once these are specified, each mass element 
goes for its own "free ride" ! 

We shall seek generalizations of the ZA by first systematically writing down a set 
of Lagrangian evolution equations for the velocity and gravity gradient for a given fluid 
element. We discuss two local approximations based on ignoring certain terms in the 
evolution equation for the Newtonian tidal tensor. One of them was introduced by 
Bertschinger & Jain (1994). They used the fact that if a quantity known as the magnetic 
part of the Weyl tensor vanishes in the Newtonian limit, the set of exact Lagrangian 
fluid equations for cold dust becomes local. (This fact was proven in general relativity by 
Matarrese, Pantano, fc Saez 1993] following earlier work of [Barnes fc Kowlingson 1959] ; part 



of the motivation for such an assumption was the statement of |Ellis 197T] that there is no 
counterpart to the magnetic part of the Weyl tensor in Newtonian theory.) Bertschinger 
& Jain then obtained the result that spindle (filamentary) collapse is favored in general as 
opposed to pancake collapse. (Pancake collapse had been thought — correctly — to be 
the generic outcome of gravitational collapse of cold dust following the work of |ZcrdovicTi 



1970[ ) Since then, it has been shown that the magnetic part of the Weyl tensor does have a 



Newtonian counterpart ( |Bertschinger fc Hamilton 1994j ; |Kofman fc Pogosyan 1995| obtained 



equivalent results but describe their conclusions slightly differently). Until now, there has 
been no quantitative determination of the magnetic Weyl term neglected by Bertschinger 
& Jain in the tidal evolution equation. With reasonable assumptions, this term may 
be negligible on super-horizon scales, leading to "silent universes" QBruni, Matarrese, fc 



Pantano 1995) 



Our second local approximation based on the tidal evolution equation is entirely new. 
It is based on dropping several more terms in addition to the Weyl tensor term. We will 
show why this is a better approximation compared to the one proposed by Bertschinger 
& Jain (1994). In fact, in tests this new approximation performs even better than the 
ZA, both in cases where exact solutions are known and where numerical solutions are 
calculated. In this paper, we concentrate on a comparison of the three local approximations 
for ellipsoids, with and without symmetries. 

To understand the main ideas underlying these local approximation methods, and how 
they differ from other approaches, it is useful to draw an analogy with gravitational lensing. 
Our use of Lagrangian fluid equations is akin to solving the optical scalar equations ( [Sachs 



1961| ), whereby one follows the two-dimensional cross-section of a congruence of light rays 



propagating through space. Our approach is similar, with light rays replaced by cold dust, 
and with the two-dimensional cross-section replaced by the three-dimensional volume of a 
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mass element. In fact, both approaches follow from the pioneering work in general relativity 
by Ehlers (1961) and Kundt & Triimper (1961). The first application of these methods to 
matter was by Hawking (1966), who pioneered the covariant fluid approach to cosmological 
perturbation theory. The formalism was championed by Ellis (1971) and eventually was 
applied to the formation of large scale structure (|Bertschinger 19951 and the references cited 
previously). As in the case of gravitational lensing, this approach can tell how a given 
(mass) element evolves but does not give its trajectory. The optical scalar equations do not 
replace the gravitational lens equation, they supplement it. Likewise, the local methods we 
advocate can supplement N-body simulations or other approximations such as Lagrangian 
perturbation theory, by providing accurate ways to follow the deformation of mass elements 
as they evolve under gravity. 

The organization of the paper is as follows. In §Q, we show how the ZA is a local 
approximation. Section |3| presents two additional local approximations based on dropping 
terms from the tidal evolution equation, and shows under what circumstances these 
approximations are exact. To compare the three different local approximations for more 
general initial conditions, in §[| we consider the motion of a homogeneous ellipsoid, in both 
cosmological (Friedmann- Robert son- Walker background) and noncosmological (vacuum) 
contexts. The Weyl tensor and other relevant terms in the tidal evolution equation are 
evaluated. In §|5| we discuss how different nonlinear approximations predict pancake versus 
spindle collapse from generic initial conditions, for which we also calculate the collapse 
times. Conclusions are presented in §|6|. The Appendix presents some results of second-order 
perturbation theory. 



In this section we review the Zel'dovich approximation starting from the Eulerian 
fluid equations in comoving coordinates. We then show that it can be regarded as a local 
approximation. 

The cosmological fluid equations for cold dust in a perturbed Robertson- Walker 



2. On the Zel'dovich Approximation 



universe with expansion scale factor a(r) are ( [Bertschinger 1995| ): 




06 





(2) 



V 2 = 4irGa 2 p6 



(3) 
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The mass density is p = p(t)(1 + 5) and v = dx/dr is the proper peculiar velocity where 
x is the comoving spatial position and r is the conformal time (hence, dr = dt/a where 
t is the proper time). We are neglecting spatial curvature so that we can use Cartesian 
coordinates where V* = V« = d/dx % for the ith spatial coordinate. 

The trajectory of a fluid element is x l (q, r) where <f is a Lagrangian coordinate labeling 
the element, conventionally chosen to be the initial position: 

x l (q,r)=q i + ^(q,r) . (4) 

Now we introduce the Lagrangian time derivative d/dr = d/dr + f- ? Vj. This time derivative 
commutes with d/dq % . Using v l = dip l /dr, we can rewrite equation (0) as 

-r^ip* + -—^ - ±irGa 2 = -VV - ^Ga 2 . (5) 
dr 2 a dr 

Each term on the left-hand side is first order in if> % . The right-hand side can be estimated 
from the Poisson equation @, but first we need the mass density. It follows in the 
Lagrangian approach by noting that pd 3 x is conserved along a fluid streamline provided d 3 x 
is computed from the mapping q — > x. If there are no displacements, q = x and p = p. The 
volume element follows from the Jacobian determinant, leading to 



dx l 



dqi 



(6) 



For small displacements the Jacobian may be expanded in a power series; 

the first-order term gives p = p(l — d^/dq 1 ) + 0(ip 2 ). Now note that 

dip 1 /dx % = (dip 1 /dq j )(dq j / dx l ) = dip 1 /dq % + 0(i/j 2 ). Therefore, using equation (Q), 

we see that the divergence of the right-hand side of equation ([5]) vanishes to first order in 

ip l . If ip l is longitudinal (i.e., has vanishing curl), then the right-hand side itself vanishes 

to first order. Displacements that grow by gravity are necessarily longitudinal in linear 

theory. The ZA consists of setting to zero the right-hand side of equation (||). (It can be 

generalized to allow for a transverse displacement; see Puchert 1993| and |Barrow fc Saich 



1993].) Under the ZA, the evolution of displacement thus obtained is used in equation 



to get the density field. The ZA is equivalent to first order Lagrangian perturbation theory 
for the trajectories x (q, r). 

With vanishing right-hand side, equation (g) is identical to the linear perturbation 
evolution equation for 5 (a fact that becomes obvious when one notes 5 = —dip l /dq l and 
d/dr = d/dr to first order in ip). This second-order ordinary differential equation in time 
has two independent solutions that we write D±(t) ( Peebles 1980|) . Taking the growing 
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solution and requiring tp l to be longitudinal, we get the solution 

mr)=D + (r) d -^ (7) 

where ip(q) is a displacement potential which is fixed by initial conditions. 

Next we will show that equations (f|) and (0) imply that the ZA displacement field is 
longitudinal in x-space (the irrotational initial conditions already imply it is irrotational in 
g-space), a first step needed before we show that the ZA is a local approximation. We have 



(V x ^ = D + (r)e ljk 



d ( dip 



dxi \dq k 




where is the usual antisymmetric Levi-Civita symbol. Now, note that the Jacobian 
matrix defined by the transformation of equations (f|) and (]7|), dx-> /dq l = 5ji+D + d 2 ip/dq^ dq l , 
is real and symmetric. By a theorem of linear algebra its inverse, dq /dx^ , is also symmetric. 
So is d 2 ip / dq k dq l and, because they commute, so is their product. Thus, in the equation 
for V x ip above, is contracted with a matrix that is symmetric in j and k, yielding 
V x $ = ( Zel'dovich fc Novikov Tggg ). 



The implication of this result is that ip is longitudinal in x-space as well as in g-space. 
The same conclusions hold for the velocity field v, since it differs from ip by only a 
time-varying factor D + /D + . As a result, under the Zel'dovich approximation we can write 

#(g(x,r),r) = D + (r) A_J > and ^(f , r ) = D+ ( T )_A_Z_i . j = , ( 8 ) 

The last equation follows from the fact that d$/dx l = dip/dq 1 (cf. eq. 0). Recall that under 
the ZA the right-hand side of equation (El) vanishes. Using equation (H), we then get 



-D + (AnGa 2 pD + Y V0 = W , (9) 



where / = ci In D + /d In a. Thus, in the Zel'dovich approximation, the velocity field is always 
(not just to first order in ip) proportional to the gravity field ( [Kofman 1991[ ). It is clear 
geometrically that this result must be correct for planar, cylindrical, or spherical flow for 
growing mode initial conditions. For plane-parallel flows, but not otherwise, the coefficient 
of proportionality of the ZA is also correct, so that the ZA is exact (until orbit-crossing) in 
one dimension. 



We are now going to present the ZA from another point of view. Similar work has 
been done by Kofman & Pogosyan (1995). Our aim is to motivate how one might improve 
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the ZA by generalizing it to a broader class of local approximations. It will become clear 
shortly exactly what we mean by local approximations. 



Let us first give a brief summary of the Lagrangian fluid equations (|Bertschinger fc 



Jain 1994|) . First of all, the gradient of the fluid velocity field is decomposed into its trace, 
traceless symmetric and antisymmetric parts, which are the expansion 9, shear and 
vorticity respectively: 

V^ = \0 5 l3+ a l3+ u l0 = u y = e yfc u; fc = (10) 

where 2u3 = V x v. Then, converting time derivatives from Eulerian to Lagrangian, equation 
(H) becomes 

^ + (1 + 5)0 = (11) 
dr 

Taking the trace of equation @ and using equations (||) and (|10D , one obtains the 
Raychaudhuri equation: 

+ °L 9 + - 9 2 + <t«<7« - 2cu 2 = -^Ga 2 pb , (12) 
dr a 3 

where u 2 = u l 0Ji. Similarly, taking the antisymmetric and traceless symmetric parts of 
equation (fj) gives respectively 

*£L + ± u * + l0u, i -o i j J = (13) 
dr a 3 J 



and 

~<h + a aij+ 3° aij + aikCjkj + UJiUJj ~ 3 Sij { aH(7kl + = ~ Eij ' ^ 
where i£y = VjVj0 — (1/3) 5ij V 2 is the gravitational tidal field. 

In keeping with the spirit of Lagrangian fluid dynamics, we would like an evolution 
equation for E^. From equations (|IJ) and @, Bertschinger & Hamilton (1994) derived 

^ + ~ Eij - V fc e kl {i H j}l + 6E i3 - + a kl E kl - 3a fe (l .^- )fe - co k {i E j)k = -AnGa 2 p a VJ . (15) 

Parentheses around a pair of subscripts indicates symmetrization, e.g., a k ^Ej)k = 
(o iEjk + a k jEik)/2. The new quantity is the Newtonian limit of the magnetic part 
of the Weyl tensor in the fluid frame. The definition and discussion of this term will be 
deferred until the next section. 

Equations ( |TT| ) to ( I5|) form a hierarchy of Lagrangian fluid equations. It is an 
incomplete set because we have not stated the evolution equation for Hij. In order to arrive 
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at a local set, we must eliminate the gradient term in equation (|T5|), either by finding an 

(i H 3)l 



approximation for — e kl uHj\i or by truncating the hierarchy in a way that eliminates our 



need to determine it. 

The ZA eliminates the need to calculate Hij by approximating the evolution of the 
gravity field — equation @ relates V0 to v. As a result, the tidal tensor in the ZA follows 
from the shear: 

AnGap 3Q H 2 
E {j = — —a* = ~2Tj-^ • (16) 

Furthermore, the divergence of the gravity field is given in the ZA by the velocity expansion 
scalar 9 instead of the density fluctuation. Thus, the ZA is equivalent to solving the local 
evolution equations 

_ + _, + _, 2 + ,. Vy -2^ = ^, (17) 

rfo-j,- d 2 n k 1 f hi 2\ 4xGap 

-^T + ^ + 3 + °" ifc(T J' + ~~ 3 °« V a + u ) = Hf ° ij ' ^ ' 

Together with equations QTT]) and (]T3|), these give a closed set of equations for the evolution 
of quantities for a single mass element with no spatial gradients. This is what we mean by 
locality. Note that we have assumed the irrotational flow initial condition and so uj^ = 
from equation (|T3| ) at all times before trajectories intersect. Equation (|18|) can be written 



also as an evolution equation for Eij by making use of equation (|16D ( |Kofman fc Pogosyan 



1995|) . But it is clear that in terms of obtaining a closed set of local equations, it is sufficient 



to stop at the level of the shear equation (18 ) 



Hence, we have shown that the ZA is a local approximation based on truncating the 
set of Lagrangian fluid equations at the shear evolution equation by setting proportional 
to <Jij and by approximating the gravitational source term in the Raychaudhuri equation. 
It is then very natural to ask whether we can go further, by using the exact Raychaudhuri 
equation and by truncating the system of equations at the tidal evolution equation with a 
different approximation from the ZA. 

There is a simple argument for why we should expect to be able to improve on the 
ZA. It is well known that the ZA gives incorrect results for spherical infall. For spherical 
infall, the velocity and gravity fields are isotropic around a point, so that cr^ = = at 
that point. Yet, the ZA overestimates the collapse time for a uniform spherical tophat. 
The reason for this is that the ZA does not obey the Poisson equation, so the right-hand 
side of equation (|17D is not exact. We can at least correct this term. We have tested this 
approximation — using equation (|T2|) in place of equation (|T7D, and using equation (pT8|) for 



the shear evolution — and found that it works poorly aside from spherical flow. Thus, we 
seek improved approximations based on a more accurate treatment of the tidal tensor. 
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3. Two Local Approximations Based on the Tidal Evolution Equation 

As remarked in the last section, the hierarchy of Lagrangian fluid equations can be 
truncated at the tidal evolution equation, provided that we approximate, or eliminate, 
the term (and possibly other terms also). If possible, we would like to find local 
approximations that retain the successes of the Zel'dovich approximation. These include 
giving the correct results in linear perturbation theory and giving the exact solution for 
plane-parallel flows. Ideally, we would also like to improve on the Zel'dovich approximation 
by giving exact results for spherical and/or cylindrical flows. We use these criteria in 
seeking improved approximations. 

Let's look at the magnetic part of the Weyl tensor more closely. The definition is given 
in Bertschinger Sz Bamilton (1994): 

H ij = - 1 -V {i H j) -2v k e k \ i E j)l (19) 

where H t satisfies: 

V x H = -16nGa 2 f ± , V • H = . (20) 
Here f± is the transverse part of the mass current, defined as follows: 

/l-i^v(^). (21) 

Using these definitions, we can rewrite equation ( |l5f) as follows: 

dE ij + a E .. + M , , = _ 4z7T Ga 2 pa ij , (22) 



where 



dr a 



My = — Vfc e kl (iHjy + OEij + Sij a kl E k i — 3a k ^Ej) k — Lo k ^Ej) k 
= —AirGa 2 pV(iVj) ^-(VjVja0) 

= -4vrGa 2 V (i f ±j) - v^VN^ + v {i V 3) V 2 . (23) 



Let us first consider plane-parallel flows, for which the ZA is exact. The velocity and 
gravity gradient tensors may be written 

V iVj = 9 diag(0, 0, 1) , = V 2 diag(0, 0, 1) , (24) 

where diagQ denotes the elements of the diagonalized tensor. Evaluating My using 
equations fl23|) , we find that the curl Hji term as well as the sum of terms proportional to 
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the tidal tensor vanish identically. The individual tidal terms do not vanish. This result 



suggests two different closure schemes for the tidal evolution equation (22). The first one 
is to discard Vfc e kl uHjy. The second is to discard the complete tensor My. If some of 
the tidal terms of M y - were retained, the resulting approximation would not be exact for 
one-dimensional flows, hence would not improve on the Zel'dovich approximation. 

The first choice, setting = in equation (|i5|), was proposed by Bertschinger & Jain 
(1994): 

^ + - Eij + 6E tJ + 6 tJ a kl E kl - 3a k {i E j)k - uo\ t E ])k = -AnGa 2 p a %3 . (25) 

We shall call this the non- magnetic approximation (NMA). Combined with equations 
(pT|)-([T^), it provides a closed set of local evolution equations. The NMA was inspired, 
in part, by the remark of Ellis (1971) that the magnetic part of the Weyl tensor has no 
Newtonian counterpart. However, it leads to unusual behavior, implying that cold dust fluid 
elements generically collapse to spindles ( [Bertschinger fc Jain 199~4|) . Also, Bertschinger 



& Hamilton (1994) were able to derive equation ( p~5|) with defined using equations 
(19|)-(12l"|) from Newton's laws in an expanding universe, as well as constraint and evolution 
equations for itself (the latter using post-Newtonian corrections), from which we now 
know that is not identically zero in the Newtonian limit, aside from some special cases 
of high symmetry. 

Thus, we are motivated to try the second approximation, setting M y - = in equation 

my- 

+ ® EtJ = -^Ga 2 pa tJ . (26) 
dr a 

Equation ( |2"6| ) and equations (PT])-(P^) form our new set of closed local equations. We shall 



call this the Local Tidal Approximation (LTA) to distinguish it from equation ([25]), the 
non-magnetic approximation. 

The LTA, like the ZA, is exact for plane-parallel flows prior to the intersection of 
orbits. What about spherically and cylindrically symmetric flows, for which the ZA is not 
exact? For the LTA, we use equations (^) to evaluate for flows that are spherically 
symmetrical around the fluid element under consideration. As long as the gravity gradient 
is finite at the origin (a condition that holds for any continuous finite-density mass 
distribution), this restriction implies that all three eigenvalues of are equal, so 

= identically (similarly = 0). Equations (p5|) and (|26|) are satisfied trivially. Thus, 
the LTA is exact for spherical mass elements. So is the NMA. 

Next we consider a non-singular fluid element on the symmetry axis of a cylindrically 
symmetric flow. By this we mean that two eigenvalues of VjVj0 are equal and the third 
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one vanishes and similarly for Vjfj. In this case we have 

V iVj = ~0diag(l,l,O) , V,V,0 = ^V 2 0diag(l,l,O) . (27) 

Using this, it is easy to show that the sum of tidal terms in the first form of equation ( f23|) 
do not vanish, while, with equations (0) and (0), the second form for leads to Mij = 0. 
Thus, the LTA is also exact for cylindrical flows, while the NMA is not. Bertschinger & 
Jain (1994) erred in saying that the NMA was exact for cylindrical flows. 

One can generalize these results to show from the second form of equation (f23|) 
that = for any flow for which (V 2 0) _1 VjVj0 equals 9~ 1 'VnVj^ and is a constant 
tensor. These conditions are equivalent to saying that the orientation and axis ratios 
of the gravitational and velocity equipotentials are constant for the mass element under 
consideration. Thus, the LTA is exact for flows with equipotentials of constant shape. 
Although this condition does not always hold, it is valid for the growing mode in the linear 
regime and it includes spherically and cylindrically symmetric flows as well as plane-parallel 
flows. Moreover, the gravitational potential contours are more nearly spherical than the 
density contours around a peak, so their shape would be expected to change relatively 
slowly with time, suggesting that the LTA may be a good approximation in general. 

In the linear regime, the LTA, NMA, and ZA all agree. It is already clear that they 
must differ in second-order perturbation theory; the Appendix presents the calculation of 
Vfe e uHjy and M^. However, it is more important to see how these various approximations 
behave as collapse is approached. We know already that generic initial conditions lead 
to collapse along one dimension (pancake) with the ZA (from [Bertschinger fc Jain 1994] ) 
while the NMA leads to collapse along two dimensions (spindle). What about the LTA? 
How accurate is the LTA for asymmetrical initial conditions? Before answering these 
questions we first examine the relative sizes of the terms in equation (^) for an overdense 
homogeneous ellipsoid in an expanding universe. 

4. Collapse of a Homogeneous Ellipsoid 

We summarize here the equations of motion for an irrotational homogeneous ellipsoid 
embedded in an expanding universe. The various interesting quantities in the tidal evolution 
equation are then calculated for the collapse of a particular ellipsoid. 

We consider an irrotational homogeneous ellipsoid with proper axis lengths Ri, R2, 
and R3 embedded in a Friedmann-Robertson- Walker background. The equations of motion 
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are ([Icke 1973fc |White fc Silk"T979l) : 



d Rj 



-2nGR, 



\pb + a iiPe ~ Pb) 



dt 2 

where t is the proper time (dt = adr) and a, is defined by 



(28) 



f°° ds 

on = R1R2R3 / = . (29) 

Jo (R4 2 + s)^J(R 1 2 + s)(R 2 2 + s)(R 3 2 + s) 

Here p e is the total density within the ellipsoid while pb is the density of the expanding 
universe surrounding the ellipsoid. They are related to the mean and perturbed densities 
used previously by pi, = p and p e = p(l + 5) . We evaluate them from the evolution of the 
axis lengths and the background expansion scale factor: 

p e R 1 R 2 R 3 = peo , PbCt 3 = Pbo , (30) 

where p eo and pbo are constants. Note that ot\ + a 2 + a 3 = 2 and we assume Oo = 1. 
Note also that since equation is second order, there are in general two independent 
modes. We choose growing mode initial conditions. For small a the first order solution or, 
equivalently, the ZA result, is 

Ri(t) = a(t)Xi (l - ^a l0 6 a^j (31) 

where the X^s give the initial axis ratios, gives the initial ellipsoid parameter, and Sq 
gives the linear amplitude of the density perturbation. We set do — 1 for overdense ellipsoids 
without loss of generality. 



Equation (pOf ) implies that the total mass, including the mass inside the ellipsoid as 
well as outside, is actually not conserved even though the mass inside the ellipsoid is. 
Hence equation (|28|) can only be an approximation to the true evolution of an initially 
homogeneous ellipsoid. In general, one expects that such an ellipsoid would cause the 
density of its immediate surroundings to deviate from the cosmic mean. Tidal fields from 
this perturbed external material should then induce departure from homogeneity in the 
ellipsoid. Based on results from an N-body simulation (S. D. M. White 1993, private 
communication), we assume that it is a good approximation to ignore departures from 
homogeneity inside and outside the ellipsoid when calculating the evolution of the axis 
ratios. 

It is noteworthy that equation ( p8| ) is exact if pb = 0, i.e., for a homogeneous ellipsoid 
in a vacuum. Later in this section we will test our approximations using the exact solution 
in this case. 
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The peculiar velocity field inside the homogeneous ellipsoid is described by: 

« = OH) * (32) 

and the gravitational potential within the ellipsoid is: 

(j) = 7iGa 2 (p e - p b ) ^2 a i x 1 ■ (33) 

i 

(The X{ are comoving coordinates and dots denote conformal time derivatives.) Quantities 
like the expansion, shear, and tidal field can be immediately read off from these expressions: 

<Tij = diag ( — - ~ TT ] , (35) 

Eij = 2irGa 2 (p e - p b ) diag fa* - . (36) 
The tensor My defined in equation fl23|) is given for the homogeneous ellipsoid by 

2\ „ 5 



Mij = lisGa p e diag 



— 1(JiA + ( Q!j 9 Cli 

3 1 3/ 1 + 5 



(37) 



Using the time evolution of Ri given by equation ( pg|) , the evolution of the various 
quantities above can be calculated. In particular, we are interested in the relative magnitude 
of various terms in the tidal evolution equations ([15]) and (|22|). We integrated equations ( p8[) 
and (|29|) numerically starting from equation (^1|) at a = 10~ 8 with axis ratios 1 : 1.25 : 1.5. 
From the axis lengths Ri and their time derivatives, using equations (|34] ) -(|3~6"1 ) we calculated 
the velocity and gravity gradient terms inside the ellipsoid. From these we then calculated 
the evolution of — e fci ^Hjy, M^, and 47iGa 2 p<Jij inside the ellipsoid. Note that in this 
test we do not integrate the tidal evolution equation itself; rather, we evaluate the terms 
in it assuming that the system evolves according to the homogeneous ellipsoid solution. 
Although, as we noted above, this solution is not exact, we are being self-consistent by 
evaluating the various tensor quantities using equations (j34l)-(j37j), which assume spatial 
homogeneity inside the ellipsoid. 

Figure 1 shows the results of this calculation. The magnitudes of — V/t e kl ^H^i and 
Mij are divided by the magnitude of 4:7iGa 2 paij, where by the magnitude of a matrix we 
mean the square root of the trace of its square. We see that — e kl uH^i and My are both 
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Fig. 1. — Evolution of the magnitudes of — Vfc e kl uHjy (dashed line) and My (solid line) 
divided by the magnitude of 47iGa 2 paij, evaluated for a homogeneous ellipsoid with initial 
axis ratios 1 : 1.25 : 1.5 embedded in an expanding universe. The magnitude of a matrix is 
defined here as the square root of sum of squares of eigenvalues. 



small compared to 4irGa 2 paij at both early and late times, but not intermediate times (near 
maximum expansion). Interestingly, the magnetic term and My have similar magnitude 
throughout the collapse process. 

We can easily understand why AirGcPpcfij is much bigger than both — Vfc e kl uHjy and 
Mij at early times using perturbation theory. The shear is first order. The last form of 
equations ( E3J ) is the best place to see that My is second order: fxi is second order because, 
to first order, fa = pVi is longitudinal (we assume irrotational initial conditions). The other 
contributions to M y - are obviously second order. From the first form of equation (23|), we 
conclude also that — e kl uHfii is second order. Expressions for these two tensors in second 
order perturbation theory are given in the Appendix. They are both nonzero in general. 



The behavior of these quantities close to the moment of pancake collapse can also 
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be estimated analytically. Suppose that the third axis collapses while the other two axes 
still have finite lengths. The term AirGa 2 paij diverges at the moment of pancake collapse 
because p diverges and so does o^-, owing to the R3/R3 term in equation (|55]). For the 
behavior of — e kl uHjy and My approaching collapse, we need to understand the behavior 
of the ctj's. 

It follows from equation ( p9|) that ol\ and a 2 vanish in the limit of vanishing R 3 for finite 
Ri and R2, because the integral is finite while the factor of R3 in front vanishes. Note also, 
by definition, the three ctj's always add up to 2. Hence a% = 2 at collapse. Moreover, it can 
be verified using equation (|29| ) that the dj's are finite at the moment of collapse, assuming 
the RiS are finite. It can then be shown using equations (|33])-([J7]) that the particular 
combination of a^'s conspires to render both Mij/A.iiGa 2 p and — e kl ( i Hj)i/A.iiGa 2 p finite. 
Hence at the moment of pancake collapse, — e kl nH^i and My are indeed much smaller 
than 47rGa 2 p<7ij. 

The fact that — e kl uH^i and My are both small compared to 4irGa 2 paij at early 
times and at the moment of pancake collapse suggests that the NMA and LTA might both 
be good approximations. However, Figure 1 shows that these terms are not negligible 
throughout the collapse process. Hence there is no guarantee that either approximation 
can reproduce the correct features of the collapse process. In particular, we do not know 
from these results whether the NMA or LTA would produce pancake collapse given the 
initial conditions we have chosen. We also do not know which approximation will be more 
accurate for generic initial conditions, although Figure 1 suggests that it may be better to 
neglect My than Hy. 

5. Pancakes Versus Spindles 

The oblate and prolate configurations are distinguished by the signature of the 
eigenvalues of and <7y. For the collapsing oblate (pancake) configuration, the eigenvalues 
of Eij have the signature (— , — , +) and those of <7y have (+, +, — ). For the collapsing 
prolate (spindle) configuration, Ey has eigenvalues with signature (—,+,+) and o-y has 
(+, — , — ). One way to see why this is true is by inspecting equations ( pop and (|36|) . For 
the pancake configuration, one can use the fact that 0.3 is close to 2 (supposing collapse 
occurs in the third direction) while a,\ and «2 almost vanish. For the spindle configuration, 
suppose that collapse occurs for the second and third direction and suppose for simplicity 
that they collapse at the same rate. Then from equation (]29f) , one can show that close to 
the spindle configuration, a 3 ~ a 2 ~ 1 and a.\ ~ 0. Using this and equations ( |35D and (|36|) , 
it is possible to obtain the signature for the eigenvalues of E^ and o"jj. Note also that it is 
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sufficient to consider only the divergent parts of £L- and cry to get the right signatures. 



Consider equation (EHp. This is the tidal evolution equation of the NMA, which ignores 
the magnetic part of the Weyl tensor. First of all, the term proportional to a/a always 
tends to decrease Ey, encouraging spherical collapse. But by the time the motion of the 
object under consideration breaks away from the expansion of the universe, this term 
becomes unimportant. Suppose now that the object is close to the pancake configuration 
with Eij having signature (— , — , +) and <7y having (+, +, — ). Then it can be seen that 
all the terms favor pancake collapse (or favor neither pancakes nor spindles) except the 
shear-tide coupling terms 5ijir kl E k i — 3a k ^Ej) k . The net sign of these two terms is such 
that the growth of E^ towards the pancake signature is suppressed. Suppose on the other 
hand that the object is close to the spindle configuration with E^ and <7y having signatures 
(— , +, +) and (+, — , — ) respectively. Then all the terms, including the shear-tide couplings, 
encourage the growth of tide towards the spindle signature. In other words, the NMA on 
the whole favors collapse toward the prolate or spindle configuration. 



Consider, on the other hand, the exact tidal evolution equation (|22|). For an object 
with a very short third axis compared to the other two, we expect «3 to be slightly less 
than but close to 2 and ai and «2 to be small and positive. Substituting this into equation 
(j37[) and looking only at the most divergent terms, one can verify that My has signature 
(—,—,+) close to the pancake configuration. Using similar arguments, it can be deduced 
that My has signature (— , +, +) close to the spindle configuration. Hence, My has the same 
signature as close to collapse, whether it be pancake or spindle; therefore it stabilizes 
collapse just like the Hubble damping term proportional to a/a. Hence ignoring My, which 
is the LTA, does not favor spindles over pancakes. This is very different from the NMA. 

Equation ( |2"B| ) tells us that My contains both — e kl uH^i and the shear-tide coupling 
terms. We can now see what is wrong with the NMA — for the spindle configuration, 
— Vfc e kl uHfii has a signature that is opposite to the shear-tide coupling terms, and it is 
large enough to reverse the spindle-enhancing effect of the latter. As a result, My as a 
whole, which includes the sum of these terms, plays no favorites. 

Numerical integration bears out this analysis. We tested the LTA and NMA by 
integrating the sets of local Lagrangian fluid equations which are obtained by ignoring the 
relevant terms in the tidal evolution equation: equations (O), (0)' (0)> an d either (^5|) or 
(f26j). The tensor equations were diagonalized along principal axes. Initial conditions were 
chosen using equations and (|34l)-(j36l) so as to correspond to the homogeneous ellipsoid 
model with initial axis ratios 1 : 1.25 : 1.5. Given the numerical solution for 9{t) and 
o"y(r), we then predicted the evolution of the homogeneous ellipsoid axis lengths by solving 
equations (|34|)-(|35"|) for Ri/Ri and numerically integrating it to get Riij). For comparison, 
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we also computed the prediction of the ZA for the axis evolution given the same initial 
conditions. We obtained the same results for the ZA by integrating the local Lagrangian 
fluid equations (|T7]) and QT%D as we did from equation (pTCf). 




a 



Fig. 2. — The evolution of axis lengths for a homogeneous ellipsoid embedded in an expanding 
universe. The initial axis ratios are 1 : 1.25 : 1.5. The "exact" solution (ignoring development 
of inhomogeneity, solid curve) is compared with the ZA (short dashed curve) and two 
local approximations: the local tidal approximation (LTA, long dashed curve) and the non- 
magnetic approximation of Bertschinger & Jain (NMA, dotted curve). 



Figure 2 compares the local approximations (ZA, LTA, NMA) for the evolution of 
the axis lengths with each other and with the solution given by integrating equations 
(p8j). Both the ZA and LTA reproduce the qualitative features of pancake collapse. As 
we have already noted, the NMA predicts collapse to a spindle instead of pancake. For 
these initial conditions, at least, the LTA is even more accurate than the ZA. The LTA 
overestimates the expansion factor at collapse by only 3%, compared with 52% for the ZA. 
The LTA appears to rectify one of the well-known problems with the ZA, namely the fact 
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that it underestimates the rapidity of collapse for non-planar perturbations. This result is 
consistent with our observation in §|3| that the LTA is exact for spherical and cylindrical 
symmetry. 



LTA=Exact 
-NMA 




.5 1 1.5 2 

a 

Fig. 3. — The evolution of the radius of a cylindrical perturbation in an expanding universe, 
corresponding to a homogeneous ellipsoid with axes R : R : oo (a cylinder). The exact 
solution (solid curve) is compared with the ZA (short dashed curve) and the NMA (dotted 
curve). The LTA is exact for this case. 



It is also useful to compare the local approximations with the exact solution for 
cylindrically symmetric perturbations. Consider a homogeneous overdense cylinder in an 
Einstein-de Sitter universe, with radius R(t). The equation of motion is given by Fillmore 
& Goldreich (1984). It can be written in a form corresponding to equation (|28|): 



d 2 R 
IF 



-2irGR 



~Pb + (fie - Pb) 



where p c is the density inside the cylinder. In fact, this is identical to equation ( 
a homogeneous ellipsoid with axis ratios R : R : oo, for which a,\ = «2 = 1, «3 



(3? 

for 
0. 
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We repeated the comparison of local approximations with the exact solution given by 
integrating equation (j5B|) . The results are shown in Figure 3. The LTA is exact, while the 
NMA underestimates the expansion factor at collapse (by 23%) and the ZA overestimates 
it (by 36%). 




a 



Fig. 4. — The evolution of axis lengths for a homogeneous ellipsoid embedded in empty 
space. The initial axis ratios are 1 : 1.25 : 1.5. The exact solution (solid curve) is compared 
with the predictions of LTA (long dashed curve) and NMA (dotted curve). 



An exact solution also exists for a homogeneous ellipsoid in a vacuum (non-expanding) 
background ( [Lin, Mestel, fc Shu 1965| ). It is easy to modify the NMA and LTA equations 
for this case, by setting a = 1 and pi, = p = 0. We did not integrate the non-cosmological 
analog of the ZA. As for Figure 2, we set the initial axis ratios to be 1 : 1.25 : 1.5, although 
in this case we set the initial velocity field to zero. Figure 4 shows the results. Once again 
we see that the LTA is rather accurate for generic initial conditions (the collapse time here 
is 1.5% too large) and leads to pancake collapse, while the NMA incorrectly predicts spindle 
collapse. 
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To compare the three local approximations (ZA, NMA, LTA) with more general initial 
conditions, we follow the notations of Bertschinger & Jain (1994) and write traceless 
symmetric tensors in terms of a magnitude and an angle: 

(Jij = \o Qij(a) , Eij = y Gpa 2 e(l + 6) {(3) . (39) 

We have introduced new scalars a < 0, e > 0, a and f3 (0 < a, (3 < n), and a one-parameter 
traceless quadrupole matrix 



Qij (a) = diag 



a + 27r\ fa — 2n\ fa 
cos | — - — , cos — - — , cos — 



(40) 



3/' V 3 J ' V3. 

With this parametrization, oblate configurations have cos a > while prolate configurations 
have cos a < 0. Of course, the shape of a perturbation can change with time. The equations 
of motion for a, e, a, and /3 for the NMA are given by Bertschinger & Jain. For the LTA, 
their equations (13) and (14) are changed to become 

d(3 3a . (a- (3\ 

One quantity of interest for general initial conditions is the expansion factor at collapse, 
i.e., the linear overdensity when a given mass element collapses. Following Bertschinger & 
Jain (1994), we parametrize the initial conditions by €q and «o, which are related to the 
values of e and a in linear theory through e = aeo and a = a$. Because initially underdense 
perturbations can collapse if the shear is sufficiently strong, we treat both initially overdense 
and underdense perturbations by specifying 5 — ±1> respectively (£ being related to 
5 in linear theory by 5 = a5o). The expansion factor at collapse, a c , is determined by 
integrating the local evolution equations for the LTA and NMA. For the ZA, it is simpler 
to use equations (f|) and (||), noting that collapse occurs when the determinant of dx l /dqj 
vanishes. With our parametrization of the initial velocity and gravity gradient tensors, it 
follows that 

a c = — — (ZA) . (43) 

S + 2e cos(a /3) K J K J 

The collapse expansion factor a c is defined to be the absolute value of the linear overdensity 

when a given mass element collapses to infinite density. For example, an overdense spherical 

perturbation collapses when a c = 1.686, while a cylindrical perturbation collapses when 

a c = 1.466 and a plane-parallel perturbation collapses when a c — 1. Although there exists 

no exact solution for arbitrary initial conditions, it is informative to compare all three 

methods. Based on our previous results we expect the LTA to be accurate to a few percent. 
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We plot contours of constant collapse time for different initial tidal parameters e and 
ao for the three local approximations in Figure 5. In each part, the left panel gives results 
for overdense perturbations while the right panel is for initially underdense perturbations. 
Figure 5b presents the same results as Figures 1 and 2 of Bertschinger & Jain (1994). 
We see that the LTA and ZA are qualitatively similar, although the ZA overestimates 
the collapse time for overdense configurations with small tide (near the center of the 
figures). According to the ZA, a c = 3 for spherical perturbations while the exact value is 
|(2/37r) 2/3 = 1.68647. . .. Both the ZA and LTA indicate more rapid collapse for initially 
oblate configurations. As noted by Bertschinger & Jain, initially prolate configurations 
collapse faster in the NMA because according to its incorrect dynamics initially oblate 
configurations must change shape before collapsing to a spindle. 

Bertschinger & Jain (1994) also noted that shear can lead to collapse of underdense 
perturbations. From Figure 5, we see that that the size of the non-collapsing region in 
parameter space (in the middle of the right-hand panels) is largest for the NMA and smallest 
for the ZA, indicating that the NMA underestimates the fraction of initial underdense 
perturbations that can collapse, while the ZA overestimates it. Using the probability 
distribution of eo and «o derived by Bertschinger & Jain for a Gaussian random field, 
we find that the probability that a randomly chosen mass element will collapse is 0.780 
for NMA, 0.888 for LTA, and 0.920 for ZA. Thus, taking the LTA as the most accurate, 
approximately 78% (= 2 x 0.888 — 1) of the underdense perturbations (and 100% of the 
overdense ones) will collapse. This estimate neglects the crossing of mass elements, which 
increases the likelihood of collapse by increasing the density. Indeed, we expect every mass 
element collapses eventually in a perturbed self-gravitating cold dust medium. 



6. Conclusion 

In this paper, we have discussed three different local approximations for gravitational 
collapse of perturbations in an expanding universe: the Zel'dovich approximation (ZA), 
the non- magnetic approximation (NMA) of Bertschinger & Jain (1994), and a new local 
tidal approximation (LTA) introduced here. Conventionally, the ZA is presented as a 
mapping of Lagrangian to Eulerian positions. However, we showed that it can also be 
regarded as a certain truncation of the set of Lagrangian fluid equations for the density, 
velocity gradient, and tide following a fluid element of cold dust. With the ZA, the gravity 
gradient is explicitly proportional to the velocity gradient, resulting in modifications to the 
Raychaudhuri and shear evolution equations. The tidal evolution equation need not be 
integrated in the ZA because the gravity field acting on a mass element is given by a simple 
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extrapolation of initial conditions. The other two approximations we discuss extend the ZA 
by integrating the exact Raychaudhuri and shear evolution equations, with approximations 
made only to the tidal evolution equation. 

All three local approximations are exact for plane-parallel perturbations. However, the 
behavior for other shapes of perturbations shows significant differences in behavior. The ZA 
is only approximate for non-plane-parallel distributions. The NMA is exact for spherical 
perturbations but not cylindrical ones. The LTA is exact for spherical and cylindrical 
perturbations and, more generally, for any growing-mode perturbations whose gravitational 
equipotential surfaces have constant shape with time. 

In order to test these approximations for non-symmetrical shapes, we compared them 
in the case of the collapse of a homogeneous ellipsoid. As expected from the results 
of Bertschinger & Jain (1994), we find that the NMA generically produces spindle-like 
singularities at collapse. The LTA, on the other hand, generically produces pancakes, 
just like the ZA. For triaxial ellipsoids, we compared numerical integrations of the local 
evolution equations with known solutions for a homogeneous ellipsoid in both cosmological 
and vacuum backgrounds. (An exact solution exists for the latter case while, in the former 
case, the homogeneous ellipsoid solution is not really exact because tides will cause the 
background, and then the ellipsoid itself, to become inhomogeneous. However, these effects 
are expected to be small.) We find that the LTA is significantly more accurate than the ZA 
(see Fig. 2). 

These results suggest we have found a promising new approximation for nonlinear 
gravitational instability. However, we have only studied the evolution of isolated irrotational 
perturbations. Caution is needed because we do not know how accurate the LTA is for more 
general initial conditions, for example, those with vorticity. Moreover, we do not know by 
how much the tide produced by other mass elements degrades the accuracy. External tides 
modulate the equipotentials surrounding a mass element; qualitatively, we expect little 
effect as long as the external tide evolves weakly or is small compared with the trace part 
of the gravity gradient. Quantitative analysis is best done using N-body simulations, which 
we leave for later work. 

The LTA has one significant limitation compared with the ZA. It tells us only the 
internal state of a given mass element (density, expansion rate, shear) and the tide on the 
element, but does not give the position of the element. However, for many purposes one 
cares more about the internal evolution of a mass element than about its position. For 
example, simple models of galaxy formation are based on spherical infall. These can be 
improved by inclusion of shear and tides (Pond fc Myers 1993| ; Eisenstein fc Loeb 1995| ). 
Our approximations could lead to even more accurate models of this sort. Also, if one 
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does need to know the positions of mass elements, then one can always supplement the 
Lagrangian fluid equations by the equation of motion for positions, perhaps using the 
Zel'dovich approximation or higher-order Lagrangian equations of motion. In principle, by 
following the velocity gradient for many mass elements, one can reconstruct the velocity 
field (up to an irrelevant overall constant), and then integrate the positions with dx/dr = v. 
An equivalent procedure was suggested by Matarrese et al. (1993). 

Perhaps the most important reason for seeking new approximations like the LTA is that 
we still lack a good understanding of the behavior of collisionless systems under nonlinear 
gravitational instability. Future work will tell whether local Lagrangian flow methods will 
provide new insights. 

We would like to thank Bhuvnesh Jain, Alan Guth, Rennan Bar-Kana, and Jim 
Frederic for helpful discussions. This work was supported by NASA grant NAG5-2816. 



Appendix: Second Order Calculation of M, 



'j 



We write the Eulerian density fluctuation field 5 = 5^' + 5^ + . . . where 5^ is treated 
as being of nth order in perturbation theory. Similar expansions are used for the velocity 
field and the scaled gravitational potential 
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For simplicity we shall assume an Einstein-de Sitter universe as in the numerical examples 
presented in this paper. In this case the perturbation series is a series in a(r). 



Peebles (1980) presents the result for 6™, which we rewrite using our variables as 
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where F 2 = F^F^ and 



F, = V,;V^ (1) 



i] — v i v ) (46) 
note that F\ = 5^. From the Euler equation (||) we get the first and second order terms 



of the peculiar velocity, 
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(47) 



where (p^ is obtained using equation ([44]) with 6™. 
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We get Mij and — e fcZ from equation (|23|). They vanish in first order; the 

second order results are 



Mr, = AixGpaa \ 5 {1) F i3 - ^V,V^ (2) + ~ [V0 (1) • V 





p. . _L -F k ,.p 



5 - {i r 3)k 



(4$ 



V fc e fc ' = il% + 4vrGpad J 3^^ - 3F fe ri F i)fc + 5 l3 [F 2 - \5<» 



(i r j)k I" Ojj 



(49) 



It can be verified that these quantities are traceless as expected using equations ( |4"5|) and 
(H3). In general, neither vanishes in second order perturbation theory. 
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Fig. 5a. — Contours of constant collapse time computed using the ZA, expressed by the 
cosmic expansion factor a c or its reciprocal, versus initial tidal field parameters. Left 
panel: initial positive density perturbations. The light (heavy) contours are spaced by 
0.1 (0.5) in a c , with the outermost contour a c = 0.4 and the central value (corresponding 
to spherical collapse) a c = 3.0. The ZA significantly overestimates the collapse time for 
low-shear perturbations. Right panel: initial negative density perturbations. The light 
(heavy) contours are spaced by 0.1 (0.5) in a" 1 , with the innermost contour a^ 1 = and 
the outermost one a" 1 = 2.3. Initial perturbations in the central region do not collapse. 
Perturbations are oblate (prolate) for eocosao > (e cos«o < 0). 
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Fig. 5b. — Same as Fig. 5a except that the NMA is used. In the left panel the innermost 
contour is a c — 1.6. In the right panel the outermost contour is a^ 1 = 1.8. The smaller 
extent of the contours for prolate configurations (e cos a < 0) reflects the fact that the 
NMA favors prolate collapse. 
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Fig. 5c. — Same as Fig. 5a except that the LTA is used. In the left panel the innermost 
contour is a c = 1.6. In the right panel the outermost contour is a^ 1 = 1.6. The LTA, like 
the ZA, favors oblate (pancake) collapse over prolate (spindle) collapse. 



